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Abstract. The Navier-Stokes equations arise naturally as a result of Ehrenfests' 
coarse-graining in phase space after a period of free-flight dynamics. This point of 
view allows for a very flexible approach to the simulation of fluid flow for high- 
Reynolds number. We construct regularisers for lattice Boltzmann computational 
models. These regularisers are based on Ehrenfests' coarse-graining idea and could 
be applied to schemes with either entropic or non-entropic quasiequilibria. We give 
a numerical scheme which gives good results for the standard test cases of the shock 
tube and the flow past a square cylinder. 
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1. Introduction 

The simulation of high- Reynolds number flow is notoriously difficult. 
In two space dimensions, a partial differential equation model for such 
flows is the Navier-Stokes equations: 
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where p, u = {u\,U2), P and E are density, velocity, pressure and 
energy density respectively. These equations model the conservation of 
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mass, momentum and energy. The number p is the coefficient of vis- 
cosity, and as this number tends to zero we recover the Euler equations 
for inviscid flow. The Reynolds number of the flow is 

Re = , 



where L is the characteristic length scale in the problem, Uqo is the 
free-stream fluid velocity, and v = p/p is the kinematic viscosity. As 
p — > 0, Re — > oo. 

Our aim is to model the flow in a physical way, so that the limit as 
the viscosity gets small is the Euler equations, but also that diffusion 
is added in a targeted, physical and controlled way. We will present a 
variant of the lattice Boltzmann method which was introduced in [7]. 
We will also set this method in the context of a more general coarse- 
graining paradigm. 



2. The lattice Boltzmann method 

Let / = /(x, v,t) be a one-particle distribution function, i.e., the 
probability of finding a particle in a volume dV around a point (x, v) , 
at a time t, in phase space is /(x, v, t)dV. Then, Boltzmann's kinetic 
transport equation is the following time evolution equation for /, 

^ + v-V/ = Q(/). (2) 

The collision integral, Q, describes the interactions of the populations 
f- 

Equation (|2|) describes the microscopic dynamics of our model. We 
will wish to recover the macroscopic dynamics, the fluid density, mo- 
mentum density and energy density. 

We do this by integrating the distribution function: 

p(x,t) := //(x,v,i)dv, 

pui(x, t) := J Uj/(x, v, t) dv, i = 1, 2, 

£(x,t) :=^J v 2 /(x,v,t)dv. 

Such functionals of the distribution are called moments. The pressure 
P is given by 

P = E- l -pu\ 
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where we have set Boltzmann's constant to 1. 

The lattice Boltzmann approach drastically simplifies this model by 
stipulating that populations can only move with a finite number of 
velocities {vi, . . . , v n }: 

9 fi 

+ Vi • V/j = Qi, i = l,...,n, (3) 

where fi is the one-particle distribution function associated with motion 
in the ith direction. 

Let m be the linear mapping which takes us from the microscopic 
variables / to the vector of macroscopic variables M: 

M := (p,pu 1 ,pu 2 ,2E). 

There are an infinite number of distribution functions which give rise to 
any particular macroscopic configuration M. Given a strictly concave 
entropy functional S(f), for any fixed M there will be a unique / which 
is the solution of the optimisation 

/^ = argmax{s(/):m(/) = M}. (4) 

We call ffij the quasiequilibrium as it is not a global equilibrium. The 
manifold of quasiequilibria, parameterised by the macroscopic moments 
M, is called the quasiequilibrium manifold. 
If the entropy is the Boltzmann entropy 

S(f) = ~ J J /log/dvdx 

the quasiequilibrium is the Maxwellian distribution 

/M(v) = ^exp(-^(v-u) 2 ). (5) 
Since m(f) = Tn(f^j) an integration rule which evaluates 

/ 5(v)/m(v) dv 

exactly for low-degree polynomials will preserve the conservation of the 
macroscopic variables M. Since the Maxwellian ([5]) is essentially Gaus- 
sian, the first candidate for this is a Gauss-Hermite-type integration 
formula. If we do this we get an approximation 

J 0(v)/(v)dv«£Wtf(vO/(vi). 
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If we write /j(x) = /(vj) then we can view the lattice Boltzmann 
equation ([3]) as a quadrature approximation in the velocity variable to 
Boltzmann's equation ([2]). For a complete treatment of this point of 
view see 1241. 



In the lattice Boltzmann community the collision of choice has be- 
come the Bhatnager-Gross-Krook [5] collision 



This is due to its simplicity and the nice intuitive interpretation that 
the dynamics relaxes towards the quasiequilibrium in a time that is 
proportional to the relaxation time r = l/u>, which models viscous 
processes. Via the Chapman-Enskog procedure it can be shown (see, 
e.g., [25J) that the associated macroscopic dynamics is the Navier- 
Stokes equations to second-order in r. 

The kinetic equation appears as an intermediate object between the 
macroscopic transport equations and the numerical LBM simulation. 
But, unfortunately, in the very intriguing limit of small viscosity and 
time step At > r, the discrete LBM model can not be a good approx- 
imation for the continuous-in-time kinetics. Nevertheless, the discrete 
model can still provide At 2 accurate approximation of the macroscopic 
transport equations [U [9]. 

We will see that it is possible to avoid the use of the kinetic equation 
as an intermediary between LBM and the hydrodynamic context. What 
we will demonstrate here (following the method presented in [151 116] ) 
is that the Navier-Stokes equations arise in a natural way via free- 
flight dynamics for a time r, followed by equilibration. The coefficient 
of viscosity will be r/2. For the remainder of the paper, the coarse- 
graining time, t, should not be confused with the relaxation time in ([6]). 
We will also show that we can approximate the Euler equations to order 
t 2 by a judicious choice of numerical scheme. A more detailed treatment 
of the construction of such numerical schemes for the approximation of 
the Navier-Stokes and Euler equations may be found in [6, 8, 9j. 



The original Ehrenfests' method [12] for introducing diffusion into a 
system was to divide phase space into cells. Then, after the dynamic 
motion of the microscopic ensemble under ([2]), which is conservative, 
an averaging occurs in the cells, giving rise to an entropy increase. 
Many other methods in statistical mechanics can be understood as a 
generalisation of this coarse-graining paradigm [13]. We will be using 




Q(f,f) 



(6) 



3. Coarse-graining 
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a modified lattice Boltzmann method to simulate the flow and we will 
describe this in the more general context of coarse-graining. 

We start from the phase flow transformation of the conservative 
dynamics: T : /(x, t) i— ► /(x,i + r). For the Ehrenfests' this was the 
flow of the Liouville equation. We mostly use the free-flight conservative 
dynamics: T : /(x, v,i) i— > /(x — vt, v,t) (this means: /(x, V, t + r) = 
/(x- vr,v,t)). 

Let r be a fixed coarse-graining time and suppose we have an initial 
quasiequilibrium distribution /q. The Ehrenfests' chain /o, /i, • • • is the 
following sequence of quasiequilibrium distributions: 

fi = /m(e T (/i_i))' * = 1,2, 

Entropy increases in the Ehrenfests' chain. By virtue of the conser- 
vative dynamics there is no entropy gain from the mechanical motion 
(from fi to @r{fi)), the gain follows from the equilibration (from T (/j) 
to fi+i). Consequently, conservative systems become dissipative. 

3.1. DETERMfNfNG MACROSCOPIC DYNAMICS 

We wish to determine the macroscopic dynamics which passes through 
the points m(/j), i = 0,1,.... In general, this will depend on the 
parameter r, so we seek an equation of the form 

Following |15}I16| we will expand this for small r in a series F(M, r) = 
Fq(M) + tF\(M) + 0(t 2 ) and match terms in powers of r to determine 
Fq and F\. In other words, for any quasiequilibrium state /q, we wish 
to have 

m(6 T (/o)) = M(r) 

to second-order in r. 

In phase space we have chosen free-flight dynamics 

f + v.V/ = 0, 

with exact solution 

G t (/ )(x,v) =/ (x-vt,v). 

Since /o is on the quasiequilibrium manifold we will replace /o with /* 
from now on. 
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The second-order expansion in time for the dynamics of the distri- 
bution / is 



e r (r) = e (f) + T 



dt 



t 2 d 2 @ t 

t=o + 2 dt 2 



t=o 



= r-rvvr+yvv(v-vr). 

Thus, to second-order, 
m(0 T (/*)) = m(0 o (/*)) - r^m(v • V/*) + ym(v • V(v • V/*)). 



Similarly, to second-order, 

dM 

M (r) = M(0) + r 



at 



+ 



- 2 a 2 M 



t=0 



2 dt 2 



t=o 



M(0) + r(F (M) + rFi(M)) + 



r 2 &F (M) 
Y dt ' 



Since M(0) = m(0 o (/* )), we have 

^2 



— rm(v • V/* 



r 



m(v- V(v- V/*)) 



= r(F (M) + rF 1 (M)) + 



t 2 dF (M) 

y 9t ■ 



Comparing the first-order terms we have 

F (M) = -m(vV/*). 
Comparing the second-order terms gives 

^i(M) + i^|^ = im(v.V(v.Vr)), 



and, upon rearrangement, we get 



1 



F\{M) = -(m(v-V(v-Vr)) 



^o(M) 
dt 



Hence, to second-order, the macroscopic equations are 

i-MH = -m(v • V/*) + J ( m(v • V(v • V/*)) - . 



at 



at a 



In what follows, to aid the flow of the presentation, we consign some 
of the calculations to an appendix. Now, let us look at the term m(v • 
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V/*) more carefully. The first component is 



mi(v- V/*) = y v V/*dv 



df* df* , 
Win h«27 — dv 



5 ■ / dv + A / «2/« dv 
7 ot 2 7 



7 5x 2 

OX\ OX2 



= V • (pu). 
Using (j24"|) . the second component is 

m 2 (v-Vf) = J uiv V/*dv 



Similarly, 

m 3 (t;-V/*) = | ^V'Vfdv 



^ ^ /n 2 

- — pnin 2 + - — [P + pu 2 
OX\ OX2 \ 



Finally, using (f25|) . 
m 4 («-V/*) = / vV V/*dv 



- / v%fdv+/- / v 2 W2 rdv 

'1 7 <7x 2 7 



(7) 



3 /n 2\ ^ 

OX 1 V / OX2 



(9) 



_d_ r o , y 

3xi - - _ - 

5 / \ 9 / \ ( 10 ) 

= 7^- (^4uiP + puiu 2 J + — ^4u 2 P + pu 2 u 2 

=i£rS ME+p) } + i-S ME+p) } 

Hence, from ([7|)- ()10p . the first-order approximation of the macro- 
scopic dynamics is (P) with \i = 0, i.e., the Euler equations. More 
specifically we have (spelling out the details for the first two compo- 
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nents), 

F 0>0 (M) = Afi - t— M 2 , (11) 

^w--^^ 3 - 11 ^^!-^^' (12) 



Mj 2 + Mf \ 


d 


Mi Mi 


M J 


dxi 


M 




d 


M 2 M 


M ) 


dxi 


M 



0,3V ; \dxtX Mo M 2 J 

9 (M 2 M 3 M 2 (M? + M$' 
+ dx~ 2 ~{ M + Mq 

We now look at the second-order correction 
' m(v • V(v • V/ )) 



2 V 9* 

Due to the computational complication of what is to follow we will only 
look at the first and second components of the above vector. This will 
give the reader a flavour of the computation. 
Let us look at the first component. Firstly, 

m (v • V(v • V/*)) = J v • V(v ■ V/*) dv 

2 2 d 2 f 
i=i j =i j 



Now, from (llip . 

9F ,o(M) 



2^ Ft AS \ A+ I 2^ Ft AS. 



dt f^ Q dMi \ dt J <9M; 



where T^jt is viewed as an operator. Now, 



9F 0fi _ 9F 0fi _ Q 9F 0fi _ d_ 8F 0fi _ d_ 

8M ~ 8M 3 ~ ' dM x ~ dxi dM 2 ~ dx 2 
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so that 

dF 0fi (M) __ d / a /_ , 2 n , d 



dt dx\ [dx 



{dx^^ + PUI ) + dx~ 2 pUlU2 } 



d { d / n\ d 



, .P + pu 2 ) + - — puiu 2 

OX2 I OX2 ^ ' OX\ 

Hence, the first component of the second-order correction is zero, which 
we would expect as this is the equation of mass conservation. 

Now we look at the second term in the second-order correction, and 
focus only on the pressure terms. 

To compute the correction terms for the Navier-Stokes equations we 
need to compute 

W ? = ' 1 ' 2 ' 3 - 

Now, we have from (fl~2l). 



0,1 



Thus, 



1 d ( M? + M$\ d MiM, 
2dxl\ ' 3 W J ~ dxi M 



dF 01 1 d 2 d 

oM 2 oxi oxi 

dF 0>1 d d 

dM\ dx\ Ul dx2 U2: 

dF 0tl d d 

DM2 dx\ 8x2 

dF 0tl 1 d 

dM 3 ~ ~2d~x[' 



Hence, 



dF ( d d 1 f 1 

mt^ = {-^ + dx- uiU2 \\-dx- pu r (13) 
dF °,i P r d 9 \ 9 ( 9 P 9 ( \ 

-dW/^ = + dx-2 U2 )dx-2 U2 \dx\ P " dx- {PUlUi) 

(14) 

dF ( <9 1 f 9 1 

-dMt F °- 2 = l&T" 2 " d^ Ul S\-d^ p - (15) 



dM 3 U,J 2 9xi I dx 
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On the other hand, from ()25[) . we have 

o 2 r 9 2 



J vwiVjf* dv = q^q x | (ui + Uj + 5ijU\jP + puimuj^ 



(^X^C^Xj J I./-' ;!/■! y 

= 3 S (UlP) + ^ (2n2P) 

(17) 

If we denote by D the difference of the pressure terms in flTTJ) and 
those in equations (fT3]) - (fT6l) then we obtain 



= 3 fi<»' p > + ^< 2 «> p > + |^ p > 

- wlsr'} " W p l + sr-W 

-^»'{^ p }- 2 fe" lP+ a^" 2P 



P) __a_„JA 

9x 2 9xi \ 9xi 



+ i_( tll P)_A Ul (Ap 

9x| 9X2 1 9X2 

9 f 9 1 9 f 9 n 



9xi 1 9x2 J 9x2 1 9xi 

_ _d_ p dui _^_ p ^± ^_ P ^1 + d p dlL2 

dx± 9xi 9x2 9x2 9xi 9x2 9x2 9xi 

= —p( —L - + —p( — + — 

9xi \9xi 9x2/ 9x2 \dx\ 9x2 

Similar calculations show that the momentum terms involving the 
derivative of terms of the form puiUj all cancel. Thus we have 

9 , x A 9 , N 9P 



J=l J 



t / _d_ p fdui _ du-A _^_ p f^l + <^± 
2 \9xi \9xi 9x2/ 9x2 \9xi 9x2 

which is the second of the Navier-Stokes equations (H|) with /i = r/2. 

Thus we have demonstrated that, in performing an Ehrenfests' step 
after free-flight we get, to the second-order, the Navier-Stokes equa- 
tions (UJ) with coefficient of viscosity r/2. This is remarkable, because 
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it does not involve any particular form for the collision integral in 
Boltzmann's equation ([2]), just free-flight and equilibration. 

3.2. Decoupling time step and viscosity 

There is of course a difficulty in simulating a Navier-Stokes flow where 
viscosity is given, with a numerical scheme in which the viscosity is 
directly proportional to the time step. The free-flight and equilibration 
scheme detailed above is such a scheme. We can write the governing 
equation in the form 

fi(*. + ViT,t + T) =//(x,t) 

= ^(x,t) + i/r(x,t), 

where /™ ir (x, t) = 2/*(x, t) — fi(x,t). Thus, after free-flight dynamics 
we move along a vector in the direction of the mirror point, / mir , which 
is the reflection of / in the quasiequilibrium manifold. With the BGK 
collision © we move some part of the way along this direction. This 
then suggests a more general numerical simulation process 

/i(x + ViT, t + r) = (1 - P)fi(x, t) + /5/r r (x, t), 

where [3 = (3{t) may be chosen to satisfy a physically relevant condition. 

A choice of (3 = 1/2 gives the Ehrenfests' step with viscosity pro- 
portional to the time step At = r. For /3 < 1/2 the viscosity is even 
bigger. Hence, in these cases the time step gives the lower boundary 
of viscosity we can realise. An important development in LBM was 
the overrelaxation step, with (3 > 1/2 (10} [T71 I26j . In this case the 
idea is that the dynamics passes through the quasiequilibrium mani- 
fold so that the next phase of free-flight would normally take us back 
through the quasiequilibrium manifold. Now this method, commonly 
called LBGK, is used for all (3 from the stability interval (3 £ [0,1]. 
For (3 — ► 1 viscosity goes to zero. One variant of LBGK is the so- 
called entropic LBM (ELBM) [18} [T9 j l2l] in which instead of a linear 
mirror reflection / i— ► / mir an entropic involution f i— ► / is used, where 
/ = (1— a)f+af*. The number a = a(f) is chosen so that the constant 
entropy condition is satisfied: S(f) = S(f). 

Both LBGK and ELBM decouple the viscosity parameter from the 
time step. There are a number of other ways in which one can achieve 
the same goal (see, e.g., [S1Q2]). We do not concern ourselves with this 
issue here but just remark that the essence is to construct a numerical 
method from the dynamics B_ T / 2 (/J^) i— > © t /2(/m)) then the first- 
order term in r is cancelled and one obtains an order r 2 approximation 
to the Euler equations. 
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After discretization of velocity space, the additional space discretiza- 
tion for LBM is not necessary in the following sense: the restriction of 
the discrete-in-time and continuous-in-space LBM chain of free-flights 
and collisions on a grid is exact, if this grid is invariant with respect to 
parallel transitions on the vectors VjT. 

Unfortunately, as we will see in Sec. H] below, there are instabilities 
in the simulation with LBGK and ELBM. This is because the free- 
flight dynamics sometimes takes us too far (to be understood in terms 
of entropy) from the quasiequilibrium manifold. In this case we apply 
a single Ehrenfests' step and return to the quasiequilibrium manifold. 
As you will see, this technique is capable of stabilising the method. In 
order to retain an order r 2 method (on average) we can only apply 
Ehrenfests' steps at a bounded number of sites. Thus we fix a tolerance 
5 which measures the entropy deviation AS = S(f*) — S(f) from its 
conditional maximum on the quasiequilibrium manifold, and then we 
choose the k (a fixed number) most distant points with AS > 5 and 
return these to quasiequilibrium. If there is less than k such points, we 
choose to return all of them. We call AS nonequilibrium entropy. 

3.3. Entropy control of non-entropic quasiequilibria 

There are several ways to define the discrete quasiequilibria /*. One of 
them is by the postulating of moment conditions: the moments m(f*) 
and their fluxes (moments of the next order, usually) should coincide 
for the discrete quasiequilibrium and for the corresponding continuous 
one (in this approach, "continuous" means "genuine"). This is the 
approach used to derive the popular polynomial quasiequilibria [25J. 
Another approach is based on an entropy condition: the discrete system 
must have its own thermodynamics and ^/-theorem, and the discrete 
quasiequilibrium should be the conditional maximum of the discrete 
entropy. 

We would like to apply Ehrenfests' stabilisation (as described at 
the end of the previous section) for all sorts of quasiequilibria. But this 
stabiliser requires the notion of entropy. In this section, we demonstrate 
how to use this entropic stabiliser for non-entropic quasiequilibria. 

Let the discrete entropy have the standard form for ideal (perfect) 
mixtures: 

S(/) = -£/*m(p|)- (18) 

After the classical work of Zeldovich [29] , this function is recognized as 
a useful instrument for the analysis of kinetic equations (especially in 
chemical kinetics 0[28]). For applications in ELBM see [20] . 
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If we define /* as the conditional entropy maximum for given 
M j = J2k m jkfk, then 

j 

where fij(M) are the Lagrange multipliers (or "potentials"). For this 
entropy and conditional equilibrium we find 

As=sun-s(f)=i^f^(j*) ( l9 ) 

if / and /* have the same moments, m{f) = m(/*). 

The right hand side of (|19p is (minus) Kullback entropy [22j . In ther- 
modynamics, the Kullback entropy belongs to the family of Massieu™ 
Planck-Kramers functions (canonical or grandcanonical potentials). 
The estimate of nonequilibrium entropy AS can be performed for both 
entropic and non-entropic quasiequilibria. Any quasiequilibrium (en- 
tropic or not) is the conditional maximum of the Kullback entropy. 
The main difference between the Kullback entropy (|19p of the form 
~J2i fi^ n (fi/ ft) an d the perfect entropy (fTHJ) is dependence of the 
denominators /* on M = m(f): f* = f* M . The perfect entropy (fTSj) is 
a free- flight invariant, and the Kullback entropy is not because of this 
dependence. 



4. Numerical Experiments 

To conclude this paper we report two numerical experiments conducted 
to demonstrate the performance of the proposed Ehrenfests' step sta- 
bilisation proposed in the previous section. The first test is a ID shock 
tube and we are interested in testing the Ehrenfests' regulariser on 
the LBGK and ELBM simulations for small (almost zero) viscosity 
(v ~ 10~ 9 ). We compare the LBGK simulation for the popular polyno- 
mial quasiequilibria [25j and for entropic quasiequilibria |20j . as well 
as ELBM simulation. In each case the scheme is supplemented by 
Ehrenfests' steps in a small number k sites with highest AS > 5. 

The second test is the 2D unsteady flow around a square-cylinder. 
The unsteady flow around a square-cylinder has been widely exper- 
imentally investigated in the literature (see, e.g., [TTJ [231 EZ])- We 
demonstrate that LBGK, with the Ehrenfests' regularisation, is capable 
of quantitively capturing the Strouhal-Reynolds relationship. The rela- 
tionship is verified up to Re = 20000 and compares well with Okajima's 
experimental data [23]. 
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4.1. Shock tube 

The ID shock tube for a compressible isothermal fluid is a standard 
benchmark test for hydrodynamic codes. We will fix the kinematic 
viscosity of the fluid at v = 10 -9 (essentially zero). Our computational 
domain will be the interval [0, 1] and we discretize this interval with 
801 uniformly spaced lattice sites. We choose the initial density ratio 
as 1 : 2 so that for x < 400 we set p = 1.0, otherwise we set p = 0.5. 

In all of our simulations we use a lattice with spacing h = 1, time 
step r = 1 and a discrete velocity set {v\, 1^2,^3} := {0, —1, 1} so that 
the model consists of static, left- and right-moving populations only. 
For a lattice site x, the neighbouring lattice sites are x + 1)2 and x + v%. 
The governing equations for LBGK are then 

f i (x + v i ,t+l) = f i (x,t) + 2p(f*{x,t)-f i (x,t)), (20) 

where the subscript % denotes population (not lattice site number) and 
fx, $2 and /3 denote the static, left- and right-moving populations, 
respectively. 

The standard polynomial quasiequilibria [25] are 

*(.-£)■ 

-(l-3n + 3u 2 ), 
6 

^(l + 3u + 3-u 2 ), 

D 

where we recall that 

P-=^2fi, pu:=J2vifi. 

i i 

For entropic quasiequilibria the entropy is S = — H, with 

H = fx log(/i/4) + f 2 log(/ 2 ) + h log(/ 3 ), 

(see, e.g., [20]). For this entropy the quasiequilibrium is available ex- 
plicitly: 

ft = ^(2-vT+8tf), 

/a =^((3n-l) + 2v / r+3^), 

/ 3 * = ~((3« + l)-2VT+3^). 



/a = 
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For our realisation of the Ehrenfests' regularisation, which is in- 
tended to keep states uniformly close to the quasiequilibrium manifold, 
we monitor nonequilibrium entropy AS at every lattice site throughout 
the simulation. If a pre-specified threshold value 5 is exceeded, then an 
Ehrenfests' step is taken at the corresponding site. Now, the governing 
LBGK equations become: 



fi(x + Vi,t + 1) 



fi(x,t) + 2[3(f*(x,t) - fi(x,t)), AS < 5, 
fi(x,t), otherwise. 

(21) 

We select the k sites with highest AS > 5 so that the Ehrenfests' steps 
are not allowed to degrade the accuracy of LBGK. 

For ELBM, entropic quasiequilibria are always employed and the 
governing equation is 

fi(x + Vi ,t + 1) = fi(x, t) + af3(f*(x, t) - fi(x, t)). (22) 

This equation differs from LBGK by the introduction of a parameter 
a which is selected to satisfy the constant entropy condition: 

S(f + a(f*-f)) = S(f). 

This is a nonlinear equation for a which we solve, using the bisection 
method, to an accuracy of 1CP 15 (see [6] for further details of the 
implementation). Supplementing ELBM with Ehrenfests' steps is the 
same as for LBGK. 

We observe that the Ehrenfests' stabilisation recipe is capable of 
subduing spurious post-shock oscillations whereas LBGK fails in this 
respect (Fig. I4.1j) . In the example we have considered a fixed tolerance 
of (k,6) = (4, 10~ 4 ). Of course, we note also that the smaller the 
tolerance 5 the more smoothing we have of the shock. Therefore we 
reiterate that it is important for Ehrenfests' steps to be employed at 
only a small proportion of the sites. 

We do not detect any advantage of using ELBM over LBGK with 
entropic quasiequilibria for this example. However, there appears to be 
some gain in employing entropic rather than polynomial quasiequilib- 
ria. We observe that the post-shock region for the unregularised LBGK 
simulations is more oscillatory when polynomial quasiequilibria are 
used. In Fig. 14.11 we have also included a panel with the simulation 
resulting from a much higher viscosity [y = 3.3333 x 10~ 2 ). Here, we 
observe no appreciable differences in the results of LBGK and ELBM. 
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Figure 1. Density and velocity profile of the 1:2 isothermal shock tube simula- 
tion after 400 time steps using (a) LBGK with polynomial quasiequilibria (|20|l 
[v = 3.3333X 10~ 2 ]; (b) LBGK with entropic quasiequilibria fl2DJ [v = 3.3333x 10" 2 ]; 
(c) ELBM (JUJ) [v = 3.3333 x 10" 2 ]; (d) LBGK with polynomial quasiequilib- 
ria (J20J) [v = 10" 9 ]; (e) LBGK with entropic quasiequilibria ([20j) [u = 10~ 9 ]; (f) 
ELBM (J52]| [i/ = 10~ 9 ]; (g) LBGK with polynomial quasiequilibria and Ehrenfests' 
steps (2TJ [v = 10 -9 , (k,S) = (4, 10~ 4 )]; (h) LBGK with entropic quasiequilibria 
and Ehrenfests' steps pi) \v = 10~ 9 , (k, 5) = (4, 10" 4 )]; (i) ELBM with Ehrenfests' 
steps [v — 10 -9 , (k, 8) = (4, 10 -4 )]. Sites where Ehrenfests' steps are employed are 
indicated by crosses. 



4.2. Flow around a square-cylinder 



Our second test is the 2D unsteady flow around a square-cylinder. The 
realisation of LBGK that we use will employ a uniform 9-speed square 
lattice with discrete velocities 



ro, 



cos((i - l)|),sin((i - 1)| 



\/2(cos(( 



. 7T 7T 

5) 2 + 4 



,sin((i - 5)| 



7T 

2 + 4 



0, 

1,2,3,4, 
5,6,7,8. 



The numbering /q, fi, ■ ■ ■ , fs are for the static, east-, north-, west-, 
south-, northeast-, northwest-, southwest- and southeast-moving pop- 
ulations, respectively. Here, we select entropic quasiequilibria by max- 
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imising the entropy functional 



s(f) = -Y,fMw) 



subject to the constraints of conservation of mass and momentum [2]: 



Here, the lattice weights, Wi, are given lattice-specific constants: Wo = 
4/9, ^1^,3,4 = 1/9 and W^e^s = 1/36. As is usual, the macroscopic 
variables are given by the expressions 



The computational set up for the flow is as follows. A square-cylinder 
of side length L, initially at rest, is emersed in a constant flow in a 
rectangular channel of length 30L and height 25L. The cylinder is place 
on the centre line in the y-direction resulting in a blockage ratio of 4%. 
The centre of the cylinder is placed at a distance 10. 5L from the inlet. 
The free-stream fluid velocity is fixed at (uqo, i>oo) = (0.05, 0) (in lattice 
units) for all simulations. 

On the north and south channel walls a free-slip boundary condition 
is imposed (see, e.g., [25J). At the inlet, the inward pointing velocities 
are replaced with their quasiequilibrium values corresponding to the 
free-stream fluid velocity. At the outlet, the inward pointing velocities 
are replaced with their associated quasiequilibrium values correspond- 
ing to the velocity and density of the penultimate row of the lattice. 
Some care should to be taken with the boundary conditions on the 
cylinder, but for more information on these the reader may consult, 



4.2.1. Strouhal-Reynolds relationship 

As a test of the Ehrenfests' regularisation, a series of simulations, all 
with characteristic length fixed at L = 20, were conducted over a 
range of Reynolds numbers The parameter pair (k, S), which control 
the Ehrenfests' steps tolerances, are fixed at (L/2, 10~ 3 ). 

We are interested in computing the Strouhal-Reynolds relationship. 
The Strouhal number St is a dimensionless measure of the vortex 
shedding frequency in the wake of one side of the cylinder: 




pu := Vi/i- 



e.g., mi. 



St 
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Figure 2. Variation of Strouhal number as a function of Reynolds. Dots are Oka- 
jima's experimental data [23] (the data has been digitally extracted from the original 
paper). Diamonds are the Ehrenfests' regularisation of LBGK and the squares are 
the ELBM simulation from [3]. 



where f w is the shedding frequency. 

For our computational set up, the vortex shedding frequency is 
computed using the following algorithmic technique. Firstly, the x- 
component of velocity is recorded during the simulation over i max = 
1250-L/uoo time steps. The monitoring points is positioned at coordi- 
nates (4L, — 2L) (assuming the origin is at the centre of the cylinder). 
Next, the dominant frequency is extracted from the final 25% of the 
signal using the discrete Fourier transform. The monitoring point is 
purposefully placed sufficiently downstream and away from the centre 
line so that only the influence of one side of the cylinder is recorded. 

The computed Strouhal-Reynolds relationship using the Ehrenfests' 
regularisation of LBGK is shown in Fig. [2J The simulation compares 
well with Okajima's data from wind tunnel and water tank experi- 
ment |23j . The present simulation extends previous LBM studies of 
this problem [3l 0] which have been able to quantitively captured the 
relationship up to Re ~ 1000. Fig. [2] also shows the ELBM simulation 
results from [3]. Furthermore, the computational domain was fixed for 
all the present computations, with the smallest value of the kinematic 
viscosity attained being y = 5x 10 -5 at Re = 20000. It is worth men- 
tioning that, for this characteristic length, LBGK exhibits numerical 
divergence at around Re = 1000. We estimate that, for the present set 
up, the computational domain would require at least ~ 10 lattice sites 
for the kinematic viscosity to be large enough for LBGK to converge 
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at Re = 20000. This is compared with ~ 10 5 sites for the present 
simulation. 
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A. Moments of the quasiequilibrium distribution 

In this appendix we calculate the moments of the distribution /* so 
as to keep the presentation in the main body of the paper clean. All 
integrals are over the whole velocity space. Firstly, we have by definition 
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Next, we have 




(23) 



Now, using the identity 




it follows by a change of variables that, for i ^ j 




Further, it follows from the identity 
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that 

|(^-^) 2 rdv = i|(v- U ) 2 rdv 

= - J (v 2 + u 2 — 2v\ui — 2v2U2)f* dv 

= ~(2P + pu 2 + pu 2 - 2pu\ - 2pu\) 
= P. 

Hence, from (|23|) . we have 

/ ViVjf* dv = SijP + puiUj. (24) 
Finally, similar calculations provide us with 

/ veViVjf* dv = (5 itj U£ + S^m + 5t,iUj)P + pu^UiUj. (25) 
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